arXiv:1505.04559v2 [cond-mat.supr-con] 13 Oct 2015 


Quasiclassical analysis of vortex lattice states in Rashba noncentrosymmetric 

superconductors 


Yuichiro Dan and Ryusuke Ikeda 

Department of Physies, Kyoto University, Kyoto 606-8502, Japan 
(Dated: October 14, 2015) 

Vortex lattice states occurring in noncentrosymmetric superconductors with a spin-orbit coupling 
of Rashba type under a magnetic field parallel to the symmetry plane are examined by assuming 
the s-wave pairing case and in an approach combining the quasiclassical theory with the Landau 
level expansion of the superconducting order parameter. The resulting field-temperature phase 
diagrams include not only a discontinuous transition but a continuous crossover between different 
vortex lattice structures, and, further, a critical end point of a structural transition line is found 
at an intermediate field and a low temperature in the present approach. It is pointed out that the 
strange field dependence of the vortex lattice structure is a consequence of that of its anisotropy 
stemming from the Rashba spin-orbit coupling, and that the critical end point is related to the 
helical phase modulation peculiar to these materials in the ideal Pauli-limited case. Furthermore, 
calculation results on the local density of states detectable in STM experiments are also presented. 
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I. INTRODUCTION 

Motivated by the recent revival of spatially modulated 
superconducting (SC) states induced by paramagnetic 
pair breakiM (PPB) M, noncentrosymmetric super¬ 
conductors [3 in nonzero magnetic fields have been stud¬ 
ied in recent years as a novel type of system with a pe¬ 
culiar modulated SC state. In noncentrosymmetric su¬ 
perconductors, the lack of spatial inversion symmetry re¬ 
sults in a splitting of the original Fermi surface into two 
sheets and makes effects of PPB anisotropic. Then, this 
anisotropic PPB effect tends to create a helical modula¬ 
tion of the phase of the SC order parameter just in a spe¬ 
cific direction 0. In particular, it is remarkable that such 
a modulated state may be realized even in a small enough 
magnetic field [2| , in contrast to the conventional Fulde- 
Ferrell-Larkin-Ovchinnikov (FFLO) states [6|, which do 
not appear unless the applied magnetic field reaches a 
high value of the order of the PPB field Hp at zero tem¬ 
perature. 

In any type-II superconductors, however, when the ap¬ 
plied magnetic field is higher than the so-called lower 
critical field iLci, the field-induced vortices enter the SC 
material, and, if the PPB-induced helical direction is per¬ 
pendicular to the applied field, the induced phase mod¬ 
ulation may be absorbed in a nontrivial manner into a 
change of the vortex lattice pattern of the SC order pa¬ 
rameter. In fact, it has been pointed out in the Ginzburg- 
Landau (GL) approach that, in the case of noncentrosym¬ 
metric superconductors with an antisymmetric spin-orbit 
coupling of Rashba type Q , any periodic phase modula¬ 
tion perpendicular to the field is gauged away in the order 
parameter solution so that the PPB-induced helical mod¬ 
ulation cannot be seen in bulk Rashba superconductors. 
However, it is unclear whether any effect of the helical 
modulation in the vortex-free limit does not occur even 
beyond the GL theory. 

In this paper, vortex lattice states in Rashba supercon¬ 


ductors which occur when the magnetic field is parallel 
to the basal plane corresponding to the symmetry plane 
for the spatial inversion are examined beyond the GL 
^proach and by combining the quasiclassical approach 
B m 5 which is wide ly exp loited in microscopic analysis 
of superconductors |1Q|-[1^ including that of multiband 
ones Bill, with the Landau level (LL) expansion of 
the order paramet er fl5l. It has been found in the previ¬ 
ous GL approach [13, [i3 fhat the structural symmetry 
of the vortex lattice in Rashba superconductors in the 
in-plane field configuration dramatically changes as the 
field increases through first-order transitions or continu¬ 
ous crossovers. This is a consequence of the enhanced role 
of the higher LLs induced by the PPB. It is difficult to 
describe such a field-dependent structural change of the 
vortex lattice in terms of the conventional method based 
on comparison in energy among a couple of assumed lat¬ 
tice structures. On the other hand, the LL expansion has 
been regarded as a convenient tool in the GL approach 
which is not applicable to lower temperatures and lower 
fields. However, the LL expansion of the order param¬ 
eter has been applied to the quasiclassical (Eilenberger- 
Larkin-Ovchinnnikov) approach to examine the diamag¬ 
netic properties [HI by incorporating an approximation 
analogous to the so-called Pesch approximation [HI . We 
have chosen to use this LL expansion method in the qua¬ 
siclassical approach to address the low-temperature 
vortex lattice states which cannot be examined in the GL 
approach 

One of the main results in the present work is the pres¬ 
ence of a critical end point of a first-order structural tran¬ 
sition of the vortex lattice in the low-temperature and 
intermediate-field regime which cannot be well described 
by the previous GL approach B- We argue that the 
presence of this critical end point is related to the heli¬ 
cal phase modulation in the vortex-free limit mentioned 
above and to the field-induced compression of the vor¬ 
tex lattice structure due to the PPB. Furthermore, as 
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FIG. 1: Fermi surfaces (FSl and FS2) under an in-plane mag¬ 
netic field. The gray arrows indicate the direction of the spin 
hxed by the spin-orbit coupling of Rashba type. Each surface 
shifts oppositely from F by ±Qo- The vector Qo is defined 
by Eq. (1^ and in Appendix A. 

an electronic measure of the structure of the vortex lat¬ 
tice at each field and temperature, we calculate the local 
density of states (LDOS) in each vortex lattice. 

The rest of this paper is organized as follows. In 
Sec. II, the electronic model examined in the present 
work and the content of our theoretical approach are 
explained. The obtained phase diagrams and the vor¬ 
tex lattice structures are shown and discussed in Sec. Ill 
together with the calculation results on the LDOS. In 
Sec. IV, our results are summarized, and the details of 
the quasiclassical analysis we have used are explained in 
Appendices. 

II. MODEL AND THEORETICAL APPROACH 

Following the previous work [T^, we start from the 
Hamiltonian with only an s-wave attractive interaction: 


H — Hsingle T Hint: 

(1) 

Hsingle ~ ^ ^ T C9k ' ^]a,yCky 


k,a,/3 


+ / d^rJ2ci{r)iisB{r)-aa,i3Cisir), 

(2) 

0.(3 



(3) 


Q 


where is the annihilation operator of an electron with 
momentum k and spin Ca{r) is that at position 

r, the cr^’s (/i = 0,1, 2, 3) are the Pauli matrices, g is the 
coupling constant, V is the volume of the system, 

— Y ^ V ^—k-\-ql2,a{~'^^2)a(3Ck-\-qj2,l3 (4) 

k,a,/3 

is the field operator of a spin-singlet 5-wave Cooper pair 
with the total momentum gr, ^single is the noninteract¬ 
ing term of the quasiparticle Hamiltonian, and Hint is 
the 5-wave pairing interaction term. Regarding the cen- 
trosymmetric term of the quasiparticle dispersion 5^, we 


assume the quasi-two-dimensional form 

Sk =-^{kl + k^) + J{1 - cos kzd), (5) 

where m is the effective mass of a quasiparticle, and d 
is the lattice constant in the c-axis direction. The anti¬ 
symmetric spin-orbit coupling (ASOC) of Rashba type is 
described by 

( 6 ) 

where = k — kzZ the two-dimensional wave vec¬ 
tor, /cp = \/2mEY^ Ey is the bare Fermi energy, z is the 
unit vector in the direction of the broken inversion sym¬ 
metry, and C is the strength of the ASOC. Throughout 
this paper, the xy plane is the basal plane on the bro¬ 
ken inversion symmetry, and J is the interplane coupling 
constant. Then, the anisotropy of the coherent lengths is 
given by 

e. V k^dJ/E^ ’ ^ ^ 

where and are the in-plane and the out-of-plane 
coherent length, respectively. The angle average over the 
Fermi surface is defined as 



for an arbitrary function h{k)^ where (j)^ = tan“^^. In 
addition, fig is the magnetic moment of the spin of a 
quasiparticle, and B is the magnetic flux density. Al¬ 
though in model (1) the orbital effect of a magnetic field 
is not incorporated, it can be readily included through 
the Peierls substitution: 

k —y k T eA, (9) 

where — e is the electronic charge, and A is the vector 
potential associated with B. 

As in the previous works SliElil , we focus on the 
case with such a realistically large ASOC that 

Max(rc, Ms 1-^1) < ICI < -E'f, (10) 

where Tc is the transition temperature at zero field. The 
smallness of the ratio Max(Tc,/rs|H|)/|C| results in sim¬ 
plifying the mean-field (BCS) Hamiltonian under which 
the Eilenberger equations are constructed. Before con¬ 
structing the mean-field quasiparticle Hamiltonian, how¬ 
ever, the quadratic term Hsingie has to be diagonal¬ 
ized. After diagonalization, we encounter a quasipar¬ 
ticle Hamiltonian consisting of two independent bands. 
In Fig.lU the resulting Fermi surfaces are sketched. 
On the other hand, this diagonalization induces pair¬ 
ing interactions between the split two bands. How¬ 
ever, these interband interaction terms are relatively of 
0((Max(Tc,/is|H|)/C)^) and hence, can be neglected. 
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Then, as explained in Appendix B, the correspond¬ 
ing transformation of the Green’s functions leads to the 
Gor’kov equations consisting only of intraband terms. As 
a result, we obtain the following Eilenberger equations: 

[2 {w„ + i(-l)“+Vs9fc ■ B} + ivp ■ n] fa 

= -2iWaAga, (11) 


and temperature, this approximated method can be used 
to determine the lattice shape to be realized over wide 
field and temperature ranges in the phase diagram. Nev¬ 
ertheless, one should keep in mind that, for the purpose 
of resolving a fine spatial structure, e.g., a single vor¬ 
tex core structure, this approximation gets less precise 
at lower fields [T5[. 

Then, Eqs. (pT]) and ([12]) are rewritten as 


[2 {cOn + i(-l)“’''Vs9fc • B} + iVF ■ n*] fa 

= -2iw*aA*ga,{12) 

ga = -^Jl + fafa < 0), (13) 

The indices a{= 1,2) specify the two split bands. Eur- 
thermore, ga, fa^ and fa are the normal and anomalous 
quasiclassical Green’s functions, t;F is the Eermi velocity 
on each ES which has the same value for both ESs up to 
the lowest order in C/^f and J/Ep (see Ref. [313 and 
also Appendix E in the present work). 


fa 9a^a-) fa 9a^a-) 9a 1 /•> ( 1 ^) 

where 

i>a = -2ma [2 {w„ + Z(-1)“+Vs9fe • -H'} • n] ^ A, 

!>„ = -2^: [2{w„ + z(-l)“+Vs9fc-Jf}+WF-n*]”^A*. 

(17) 

The order parameter A can be expanded in terms of LLs: 

A{r) = '^dNtpN{r), (18) 

N 


Wa =(-l)“iexp(i(-l)“())fe) 

={-l)°'i\\i±\~^{kx+i{-l)°'ky) (14) 

is a pairing function associated with the two bands oc¬ 
curring after the diagonalization. 


where 


X — mu) (19) 


IT — —iV T 2eA., 

n*= iV + 2eA, (15) 

and ujn{> 0) is the fermion Matsubara frequency. This 
set of equations is equivalent to that used in previous 
works [T^, lU, except for the inclusion of the Zeeman 
effect. Note that the two bands split by the ASOG are 
coupled with each other only through the shared order 
parameter A due to the condition C ^ Tc. 

Next, to solve Eqs. m and (|12p . we assume the type- 
II limit hereafter so that B = where H is the applied 
field along the y axis. In addition, following Adachi et 
al. [IBI, the Landau level (LL) expansion of A is used in 
the quasiclassical approach. This is because the conven¬ 
tional treatment based on comparison in the free energy 
among a couple of assumed structures is not fruitful in 
the present issue where field-dependent and continuous 
changes of the vortex lattice structure are expected to oc¬ 
cur [ll, [l3- Nevertheless, it is difficult to find an exact 
solution of Eqs. m and m based on the LL ex pan sion 
method, and hence, we adopt an approximation [l5| un¬ 
derestimating spatial variations of |Ap to be included in 
the normal Green’s function ga corresponding to an ana¬ 
log of the Pesch approximation M- The result following 
from this approximation, dubbed the “approximate solu¬ 
tion” in Ref. [lH, was argued there to be valid not only 
near the Hc 2 line but also in the low-field region as long 
as thermodynamic quantities are considered [2^. Since 
the central part of our present work is to find the vortex 
lattice structure with the lowest free energy at each field 


is the Ath LL {N > 0) when the Landau gauge A = 
—Hxz is chosen. Here, vh = lly/2eH is the magnetic 
length which characterizes the spacing between two vor¬ 
tices. 

The wave vector of the helical phase Q is nonvanishing 
as far as dN is finite. Throughout the present work, Q 
is fixed to 2SNQo^ where 

Qo = l^x (20) 

vf 

is the shift of the Eermi surfaces (see Eig.[T]and Appendix 
A), and the deviation of the true Q from 2SNQq is as¬ 
sumed to be compensated by incorporating as many LLs 
as possible. It is originally known that the identification 
Q = 2SNQci is safely valid near Hc 2 {T) at high tempera¬ 
tures [2 ,[i 3 where the higher gradients may be neglected. 
In the GL approach in Ref. [l7| , the validity of this iden¬ 
tification has been tested in the simplest 5-wave pairing 
case by comparing with the exact result obtained by de¬ 
termining the Q value minimizing the free energy at each 
magnetic field, and the simplified treatment based on the 
identification Q = 25NQq has been shown not to affect 
the phase diagram even quantitatively (see Eig. 3 and its 
related discussions in Ref. [13). 

The function Tat is expressed by the Ath Hermite 
polynomial Aat as follows: 


'^n{x) 


V2^A!7ri/4 ■ 


( 21 ) 
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Parameters v and A, which represent the lattice shape, 
are defined as in Fig.[2l Due to this expansion, the func¬ 
tion can be described as the linear algebraic expres¬ 
sion: 


= -2iWa ^ IpMMM^dN, (22) 

M,N 

and can be calculated from the relation, 

^a{k) = ^a{-kr. (23) 

which can be proved with the symmetry relations pre¬ 
sented in Appendix D. The matrix is defined as 

follows: 

pOC 

MIin= / 

Jo 


which is remarkably independent of the lattice shape, 
where 


s = 


+ ij 


(25) 


min(M,Ar) 

jjmn{z) = 


^/mW\ 


1=0 


(M-/)!(A-/)!/! 


(26) 


The derivation of Eq. (j22j) is shown in Appendix G. 

Since the shape of the vortex lattice to be realized is 
determined through minimization of the free energy per 
unit volume, we need an expression of the free energy 
represented by the quasiclassical Green’s functions. Ac¬ 
cording to the theory of Eilenberger [sl , the free energy 
F is calculated through the variational principle in the 
form 


= Nj 


d-^r 


|A|2lnF + 27rT ^ ^ 

u;n>0a=l,2 


1 + {-IY5N / ( |A|2 iWaMa + iA*</„ 


UJji 


9a - I 


FS 


(27) 


where T is the temperature. 


on the derivation of Eq. (|27)) are given in Appendix E. 


N = 


N 1 FN 2 

2 


(28) 


N 1 +N 2 



(29) 


and Na is the normal DOS on the ath ES. More details 


To make the optimization feasible, the LLs are trans¬ 
formed into the linear combinations of the LLs by diag¬ 
onalizing the quadratic term, F 2 , of the free energy with 
respect to the order parameter A. Using the parameter 
integral A~^ = dp exp(—pA), the quadratic term F 2 

is easily obtained from Eq. (EZl) and becomes 


V 


rj-i POO 

N du ^M,N In — + / dpf{p) {Sm,n 
M,N L c Jo 


dN 


(30) 


which coincides with the expression obtained in the pre¬ 
vious study Cl , where 


f{p) = 


2'kT 

sinh(27rT p) 


(31) 


The matrix to be diagonalized is the expression between 
the square brackets in Eq. (l3Q]) . and the resulting modes 
are the linear combinations of the LLs. If the modes re¬ 


sulting from the diagonalization are separated in energy 
from one another, we only have to select just the mode 
with the lowest energy to obtain the vortex lattice struc¬ 
ture in equilibrium, because an energy difference between 
lattice structures is usually much smaller [^. Thus, we 
have three variational parameters: A, and the ampli¬ 

tude of the relevant mode. 

After the diagonalization, the Hc 2 line is determined. 
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as usual, as the line in the H-T phase diagram on which 
the eigenvalue of the lowest energy mode changes its sign 
upon cooling. As already noted elsewhere [1^, the 
transition at the Hc 2 {T) line defined in the mean-field 
approximation is, in contrast to that in the centrosym- 
metric case [l|, of second order irrespective of the tem¬ 
perature and the strength of the PPB effect. Indeed, we 
have confirmed that the quartic term in F with respect 
to A is positive for the mode with the lowest eigenvalue 
even in the high-held and low-temperature region when 
SN = 0, where the suppression of the Zeeman effect due 
to the ASOC is the weakest. 

Furthermore, as a physical quantity testable in STM 
experiments and rehecting the vortex lattice structures, 
we have considered the LDOS. To obtain this quantity, 
the analytic continuation, iujn E if] is performed, 
where 77 is an inhnitesimal and positive. In the present 
formalism, this is equivalent to the replacement iuJn 
E if] in which leads instantly to the retarded 

quasiclassical Green’s function g^. Then, we have 

N{r-,E) = -7V^(1 + (-l)“5iV)Re 

a 

(32) 

as the LDOS (see Appendix H). In this case, note that 
the relation between the retarded versions of and 
is given not by Eq. (|23]) but by 

^a{k,E) = ^a{-k,-Ey. (33) 

III. RESULTS 

In this section, our calculation results on the phase di¬ 
agram are shown and explained. In all of our calculation 
results presented in this paper, we have commonly used 
the parameter values J/Ep = 0 . 1 , Hy^/Hp{cc gg) = 2 . 0 , 
and d = Tr/fep, which are the same values as those in 
Fig. 3 in Ref. [l^, where (= O.560o/27r^o) 

(= 1.25Tc//is) are the orbital pair breaking field in 2D 





FIG. 2: Definition of the parameters u and A, which represents 
the shape of the vortex lattice. Here, ai and a 2 are principal 
lattice vectors of the lattice. 



T/R 

FIG. 3: (Color online) Dependence of the Hc 2 curve on the 
number of LLs incorporated in the calculation in the 6N = 0 
case, where rimax is the index of the highest LL incorporated. 


systems and the paramagnetic pair breaking field at zero 
temperature, respectively, 0o = 7r/e is the flux quantum, 
and ^0 = vfI2'kTc is the coherent length in the directions 
parallel to the basal plane. 


A. (5A = 0 case 

In this subsection, the resulting phase diagram in the 
limiting case with dN = 0 is explained. According to 
the inequality m, this case corresponds to the limit of 
a large bandwidth. 

First, the number of the LLs to be incorporated in 
our calculation should be determined. As more LLs are 
included, the resulting Hc 2 value at each temperature 
becomes higher. In Fig.O such an example of the de¬ 
pendence of Hc 2 {T) on the number of the incorporated 
LLs is presented, where rimax is the index of the high¬ 
est LL incorporated. Ideally, the saturation of Hc 2 value 
should be reached by a finite value of nmax- Based on the 
^max dependence of the Hc 2 {T) curve obtained in Fig.[3l 
we have kept just the lowest eight LLs to determine the 
vortex lattice structure, as in Ref. [I 3 . 

Next, the details of the mode splittings resulting from 
the diagonalization in F 2 are explained. In Fig. HI the 
field dependencies of the eigenvalues of the diagonalized 
modes are shown at a low temperature. In this 6N = 0 
case, the lowest two modes are found to be nearly degen¬ 
erate for H > Thus, the free energies resulting 

from the two modes have been calculated individually 
and compared with each other to determine the vortex 
lattice structure in equilibrium. 

The resulting phase diagram is shown in Fig.O in 
which there are three phases (I-III) separated by first- 
order structural transitions (FOSTs). The structure in 
phases I and II may be regarded as a stretched triangu¬ 
lar lattice. On the other hand, a crossover from a one- 
dimensional-like structure, in which a vortex layer and 
a nodal line are alternating, to a honeycomb-like vortex 
lattice occurs in phase III. In the former (low-field) struc¬ 
ture of phase III, the alternation occurs along the x axis. 
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H/H^Pb 


FIG. 4: (Color online) Field dependence of eigenvalues of 
the modes obtained by diagonalizing F2 at T = O.lTc when 
SN = 0. The modes with the lowest two eigenvalues (the red 
and blue lines) are nearly degenerate for H > OAH^^ with 
each other and cross at = 0A8H^^. 


namely, the direction of the shift of the Fermi surfaces. 

The most remarkable character seen commonly in 
these vortex states is the lattice compression parallel to 
the X axis occurring with the field increasing. This is 
a consequence of the shift of Fermi surfaces caused by 
the in-plane applied field due to the interplay between 
the ASOC and the PPB effect. Without PPB, no such 
field-induced anisotropy arises. To avoid any confusion, 
the length scales in the x and z directions are measured 
in units of and respectively, hereafter. 

Reflecting the above-mentioned shift of the Fermi sur¬ 
faces, the system favors the periodicity proportional to 
l/iQ^rn) oc along the x axis in real space. On 

the other hand, since, due to the flux quantization, the 
area of the unit cell of the vortex lattice (in the above- 
mentioned units) is kept constant, the lattice spacing par¬ 
allel to the z axis is expanded with the field increasing. 
In this manner, the field-induced lattice compression in 
the X direction is explained (see also Fig. [6]). 

This field-induced compression parallel to the x axis 
tends to induce FOSTs between different vortex lattice 
symmetries in different ways. In general, a compression 
in one direction merely enhances the anisotropy of the 
lattice, and a superfluous compression accompanied by 
no change of the lattice symmetry would lead to some 
energy cost. Then, a FOST to a more isotropic lattice 
state may occur. When a couple of lattice symmetries are 
competitive in energy to each other, however, it is possi¬ 
ble for a FOST to occur between the two states without 
releasing the anisotropy. In the present bN = 0 case, the 
FOST of the latter type seems to be realized between the 
I and II phases in the low-field regime where the vortex 
lattice solution is formed in terms of only the LLs with 
even indices. On the other hand, in higher fields, the 
vortex lattice consisting only of the LLs with odd indices 
is favored because of enhanced roles of the PPB in higher 
fields. The above-mentioned release of the anisotropy is 
realized through the FOST between I and III phases to¬ 
gether with this switching in the description of the order 



T/Tc 
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FIG. 5: (Golor online) Resulting phase diagram and vortex 
lattice states appearing when 8N — 0. The first-order struc¬ 
tural transition (FOST) points are numerically determined, 
and the line connecting between them is a guide to the eye. 
Phases I and II are stretched triangular lattice phases. Phase 
III is a modulated triangular lattice phase. (B) shows that 
the modulation is along the direction of the shift of the Fermi 
surfaces. The figures (G), (E), and (F) indicate that the lat¬ 
tice is compressed along the x axis with increasing field. In 
(A), the spatial modulation in the region surrounded by the 
broken white circle does not imply the presence of an addi¬ 
tional vortex there. Here, Ao is the magnitude of the order 
parameter at each temperature in the absence of magnetic 
helds. 


parameter from the even to odd LLs. 

On the other hand, at a glance, one might wonder why 
no FOST occurs between the structures (A) and (B) if 
noting the appearance in (A) of an additional modulation 
of the order parameter amplitude indicated by the bro¬ 
ken white circle. However, this modulation suggesting a 
node of the order parameter is not accompanied by any 
nonvanishing winding number and thus, is not a genuine 
vortex but just a modulated structure with a low but non¬ 
vanishing amplitude of the order parameter. Hence, the 
structure change between (A) and (B) can occur grad- 
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FIG. 6: (Color online) Field dependence of z/, which is the 
lattice spacing along the x axis (see Fig.[2|), of the lattice in 
phase I of Fig.[5l Note that ly is almost proportional to 1/\/H 
for H > O.SFT^rb at T = 0.3Tc. 

ually and continuously to compensate the anisotropy of 
the vortex lattice as the field increases, which leads to the 
continuous crossover between them rather than a FOST. 

As mentioned above, the appearance of the vortex lat¬ 
tice consisting only of odd LLs in higher fields stems from 
the PPB effect, and consequently, the resulting vortex 
lattices in higher fields are mostly occupied by the spa¬ 
tial regions in which the order parameter amplitude |A| 
nearly vanishes. To correctly describe such vortex lat¬ 
tices with PPB-induced additional modulati ons on the 
length scales of the magnetic length vh = y/0o/(27ri7), 
the nonlocality needs to be taken into account properly in 
the terms distinguishing different lattice structures in the 
free energy. In the previous works based on the GL free 
energy kept up to the quartic order in the order param¬ 
eter A, the quartic term has been assumed in a spatially 
local form. In describing details of the lattice structure, 
this local form is insufficient particularly in higher fields 
where the PPB effect is not negligible. In fact, the result¬ 
ing structures in the phase III are different from those in 
the previous GL approach. On the other hand, the field- 
induced transition from a triangular structure to another 
one, namely, from I to II in Fig. [5l is qualitatively similar 
to the previous one [mil. 


B. SN = 0.1 case 

Next, we turn to a more realistic case with a nonvan¬ 
ishing SN or a finite bandwidth. The choice of the value 
SN = 0.1 seems to be reasonable if one images the mate¬ 
rials including CePtsSi a as the corresponding model 
systems. 

As in the SN = 0 case, the phase diagram is examined 
by including the lowest eight LLs (see Fig.f^. In this 
case, the even and odd LLs are mixed in every mode 
resulting from the diagonalization. Furthermore, as seen 
in Fig. [SI there is no competition between the modes, 
and we have a well-defined mode with the lowest energy 
eigenvalue. Thus, we only have to focus on this mode to 



T/T, 

FIG. 7: (Golor online) Dependence of the Hc 2 curve on the 
number of LLs incorporated in calculation when SN = 0.1. 
As in the SN = 0 case, we have assumed Umax = 7. 





FIG. 8: (Golor online) Field dependence of the eigenvalues of 
the modes obtained by diagonalizing Fb at T = O.lTc when 
SN = 0.1. The mode with the lowest eigenvalue (red line) is 
dominant at any field. 


determine the vortex lattice structure at each field and 
temperature based on the free energy ([27|). 

Figure|9] shows the resultant phase diagram and vortex 
states. There, the phase I and II are stretched triangular 
lattice phases, while a rectangular lattice is stable in the 
phase III. The phase IV is characterized by a modulated 
triangular lattice structure. The most remarkable differ¬ 
ence of this phase diagram from Fig. [5] in the SN = 0 case 
is the emergence of a critical end point of the FOST line 
between phases I and 11. 

First, to elucidate the effect of finite SN^ as in the 
former SN = 0 case, the relation between the lattice 
spacing along the x axis and the applied field has been 
plotted. It is remarkable in Fig. [10] that, although the 
spacing in the x direction, broadly speaking, shrinks with 
the field increasing reflecting the compression induced by 
the ASOC, an upturn appears in its field dependence. 
This seems to result from the appearance of another pe¬ 
riodicity caused by the emergence of the helical phase in 
the vortex-free case [Ij]. The wave vector Q of the he¬ 
lical phase is known to be proportional to SN (see the 
description below Eq. ([T9]) or Refs. [I, [l3): and thus, it 
is natural to expect that the effect of Q becomes larger 
as SN increases. In general, the magnitude of Q is not 
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FIG. 9: (Color online) Phase diagram and vortex lattice 
states when SN = 0.1. The first order structural transition 
(FOST) points are numerically determined, and the line con¬ 
necting between them is a guide to the eye. Phases I and II are 
stretched triangular lattice phases, while phase III is a rect¬ 
angular lattice one. Furthermore, phase IV is a modulated 
triangular lattice phase. Figure (A) shows that the modu¬ 
lation develops along the direction of the shift of the Fermi 
surfaces with the field increasing. (D), (E) and (F) can be 
continuously transformed to one another circumventing the 
critical end point (CEP). 


commensurate with that of Qq^ and hence, the role of Q 
can interfere with that of Qo, which is thought to lead to 
the upturn. 

Furthermore, the emergence of the critical end point 
of the low-field FOST seems to be closely related to the 
upturn in Fig.[T 0 l It is found in our calculation that, in 
the parameter (z^. A) space, where u and A are defined 
in Fig. [2] as parameters characterizing the vortex lattice 
unit cell, phases I and II correspond to two neighbor¬ 
ing valleys to each other. In the previous 6N = 0 case, 
these two valleys move to the same direction in the pa¬ 
rameter space as the field, and hence Qo, increases, and 
consequently, they do not merge with each other. On the 
other hand, in the present case, the effect of the nonzero 
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FIG. 10: (Color online) Field dependence of u, the lattice 
spacing parallel to the x axis, of the states in phase II of 
Fig.O The upturn of the width for OAH^^ < H < 
means that the field-induced compression parallel to the x 
axis is weakened in higher fields. As for the definition of ly, 
see Fig. [21 


Qo is weakened by the presence of the finite Q in par¬ 
ticular at lower temperatures, and hence, the structure 
in phase II, which is more strongly compressed in the x 
direction than that in the phase I, starts to return to 
a more stretched structure at a (5A^-dependent value of 
the applied field. Therefore, the above-mentioned two 
valleys tend to merge with each other, resulting in the 
disappearance of the low-field FOST between the I and 

II phases and thus in the critical end point. In fact, the 
change of the lattice structure in Fig. [9] from (E) to (D) 
and then to (F) can be naturally understood as being 
due to the field-induced compression in the x direction 
and stretch in the 2 : direction. 

In the higher-held region where the anisotropic trian¬ 
gular lattice is destabilized, the resulting structure (C) 
has the rectangular symmetry. Interestingly, this phase 

III with the rectangular symmetry is wide, and, with no 
transition, the nodelike region with extremely small |A| 
becomes wider as the held grows. This increase of the 
spatial modulation of |A| in phase III is a consequence 
of the roles of the odd LLs due to the enhancement of 
the PPB ehect in the higher-held region. Furthermore, 
at the high-held end, we have the narrow phase IV with 
highly anisotropic and modulating structures. 


C. Calculation of LDOS 

As available results for comparison with real experi¬ 
ments to be performed in future, the LDOS of vortex lat¬ 
tices have been examined, and their examples are shown 
in Fig.fTTl The smearing factor 77 is hxed at the value 
where the 77 dependence of the spectrum of the LDOS is 
moderate. 

In the low-held region, there is a double peak structure 
with a narrow splitting around E" = 0 in the vortex core, 
and as the held increases, the splitting of the peaks grows 
wider. This splitting seems to stem from the Zeeman 
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effect [ill , because its width is nearly equal to the double 
oi fjisH (= 0.4 • • 27rTc in the present cases). 

Although the peaks are split due to the Zeeman effect, 
the spatial dependence of the LDOS at ^ = 0 reflects 
that of |A| directly. Thus, observation of these peaks 
around the vortex cores in STM experiments would lead 
to the verification of the compressing effect due to the 
finite Qq (i.e., the shift of the split Fermi surfaces induced 
by the interplay between the PPB effect and the ASOC). 

We should comment, however, on the smearing fac¬ 
tor T]. In the present approach, the value of is 

of 0 ( 10 “^) and is much larger than that in the meth¬ 
ods used in such papers as Ref. 0 and Ref. 0 , where 
77/27rTc = 0(10“^). Therefore, the results here are highly 
smeared, and the detailed information on the electronic 
structure may be lost. Nevertheless, we believe that the 
essential structure is captured because there is a good 
correspondence between the vortex core and the peak of 
the LDOS. 

We may be able to overcome this difficulty within the 
present framework, where the Eilenberger equation and 
the Landau level expansion are combined, by calculating 
the “full solution,” which is introduced in Ref. M- As 
is mentioned in Sec. im the validity of the approximate 
solution is not ensured in computing quantities related 
to the fine spatial structure of the system, which might 
result in the large 77 . We can get over this point by ex¬ 
amining the full solution, where the Fourier transform of 
the normal quasiclassical Green’s function is employed 
instead of the approximation analogous to that used by 
Pesch 0 - Development of the method to calculate the 
LDOS in this direction may be done in future works. 


IV. SUMMARY 

In this paper, possible vortex lattices in Rashba non- 
centrosymmetric superconductors under magnetic fields 
parallel to the basal plane have been studied based on 
the quasiclassical approach, and the obtained results have 
been compared with those in the mevious GL approach 
00 neglecting the nonlocality [l| in the quartic term 
of the GL free energy. 

We have found that the overall field dependence of 
the vortex lattice structure in the realistic SN 7 ^ 0 case 
remains unchanged even in the quasiclassical approach: 
The lattice structure is hexagonal in lower fields, while it 
is rectangular in higher fields. 

However, we have also noticed that the details of the 
field dependence of the lattice structure are significantly 
changed. First of all, the lattice structure at the high- 
held end is signihcantly changed compared with the cor¬ 
responding GL result because of our proper treatment of 
the nonlocality of the quartic term in the free energy with 
respect to the superconducting order parameter. The 
need for such treatment comes from the strong PPB ef¬ 
fect. Due to this effect, higher LLs, the effect of which is 
known to make the lattice structure complicated in the 




-4 -2 0 2 4 




FIG. 11: (Color online) Panels (a) and (b) are the images of 
the local density of states (LDOS) at the excitation energy 
F = 0 of the states (F) and (D), respectively, in Fig.[9l The 
smearing factor 77 here is set at 0.1 and 0.15, respectively, in 
units of 27tTc. Panels (c) and (d) are the graphs of the energy 
dependence of the LDOS at the vortex center of the same 
state as in (a) and (b), respectively, with different smearing 
factors 77 . In (c), the spectrum is stable around 77/27rTc = 0.1 
and has slightly split peaks. In (d), the spectrum is stable 
around 77/27rTc = 0.15 and has widely split peaks. 


context of centrosymmetric superconductors [l^, play a 
more important role in the higher-field region. Mean¬ 
while, their spatial variation is more intense than that 
of the lowest LL. Therefore, the nonlocality has to be 
dealt with appropriately in such region. Second, we have 
found a critical end point of a first-order structural tran¬ 
sition line at a low temperature and an intermediate field, 
where the previous GL approach 00 did not give any 
reliable result. Furthermore, its appearance has been ar¬ 
gued to be a reflection of the helical phase modulation 
which can be directly seen only in the vortex free limit. 

Moreover, in the present work, we have been able to 
clarify that the origin of the complex field-dependent 
structural changes of the vortex lattice consists in the 
anisotropic compression of the lattice occurring as a con¬ 
sequence of the relative shift of the two Fermi surfaces 
due to the interplay between the PPB and the lack of 
the inversion symmetry. In addition, we have also ex- 
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amined the LDOS in such vortex lattice structures. We 
hope that, through some STM experiments, the strange 
field dependencies of the vortex lattice structure would 
be verified by measuring the LDOS in Rashba supercon¬ 
ductors. 
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Appendix A: Diagonalization of T^singie 

For convenience, we define here again the noninteract¬ 
ing part of the Hamiltonian ^single and its concomitant 
quantities: 

"^single ~ ^ ^ Cdk ' ^]q;,/3CA;/3 

k,a,/3 

+ jSr^c^^{r)^isB{r) ■ (Toc,pcp{r), (Al) 

where Cka is the annihilation operator of an electron with 
momentum k and spin CQ,(r) is the counter¬ 

part at the position r in the real space, and the cr^’s 
(/i = 0,1,2,3) are the Pauli matrices. As to the cen- 
trosymmetric part of the quasiparticle dispersion 5^, the 
quasi two-dimensional form is assumed: 

ek =-^{kl + kl) + J(l-coskzd), (A2) 

where m is the effective mass of a quasiparticle, and d is 
the lattice constant in the c-axis direction; 

(A3) 

is the {g) vector expressing the antisymmetric spin-orbit 
coupling (ASOC) of Rashba type, k± = k — kzZ is the 
two-dimensional wave vector, = \/2mEY^ Ep is the 
bare Fermi energy, z is the unit vector in the direction of 
the broken inversion symmetry, and is the strength of 
the ASOC. Throughout this paper, the xy plane is the 
basal plane for the broken inversion symmetry, and J is 
the interplane coupling constant. In addition, jig is the 
magnetic moment of the spin of a quasiparticle, and B 
is the magnetic flux density. 

In the absence of magnetic fields, ^single can be diag¬ 
onalized by using the matrix 

TT(U\ CTO+i(cos + sin 0^(72) 

U {k) = -- (A4) 

(0fc = tan“^^), namely, by introducing the new field 
operator 

Cafe — Uaaij^^^otk' 


Then, ^single with B = 0 becomes 

^single ~ ^ £2k^2k^‘^k} 

fe 

where 

Sak=Sk + {-l)°-+^C\gk\. (A7) 

In the presence of a homogeneous magnetic field iT, 
the existence of the Zeeman energy term does not per¬ 
mit the precise diagonalization by using U{k). Never¬ 
theless, if the temperature T and the magnitude of the 
Zeeman energy ggH are sufficiently small compared to 
the strength of the ASOC C, the interband mixing caused 
by the Zeeman term can be quantitatively neglected, so 
that ^single simply becomes 

^single — X {(eifc sQk ’ 
k 

+(c2fe - ys9k • iT)4feC2fe|, (A8) 

in which the momentum Qq giving the shift between the 

two Fermi surfaces is given by 

Qo = (A9) 

VF 

because 

g.Qk (AlO) 

VF 

is satisfied near the Fermi surface, where iJp is the Fermi 
velocity. 


Appendix B: Derivation of the Eilenberger Equation 

Throughout the present paper, we use the mean field 
approximation 

E ^ ^single — — + ~ l^qP 

q ^ q 

(Bl) 

for the Hamiltonian 'H, where g{> 0) is the coupling con¬ 
stant, V is the volume of the system, and 

^q ^ ^ ^ ^—k-\-q/2,a{ '^^2)al3^k-\-q/2,g 

k,a,/3 

is the field operator of a spin-singlet s-wave Cooper pair 
with the total momentum q. The component Aq with 
the momentum q of the order parameter A(r) is given 
by 

= (B3) 
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where A(r) = ^ ^Q' Here, (A) is the grand 

canonical ensemble average of an arbitrary operator X 
under the Hamiltonian H. 

As usual, the Gor’kov Green’s functions are defined as 

Ga/3{ri,r2-,Ti -T2) = - (TrCa(ri,T2)c/3(r2,T2))g^ , 

Ga/3{ri,r2-,Ti -T2) = - (IVCa(ri,ri)C/3(r2,T2))gq , 

Fal 3 {ri,r 2 ]Ti -T2) = - {TTCa{ri,Ti)cp{r2,T2))^g , 

Fap{ri,r 2 -,Ti - T 2 ) = - {TrCa{ri,Ti)cp{r 2 ,T 2 )) , 

(B4) 

_I 

-driGap{ri,r2;Ti - T 2 ) =<5^(ri - r2)<5(ri - T2)5a0 

+ +eA(ri)) + UsO- • B(ri)]„^G^/3(ri, r2; ti - T2) -E Aq;^ (ri)F^/ 3 (ri,r 2 ;Ti - r 2 ), 

7 7 

-5TiGa,3(ri,r2;ri - r2) =S^{ri - r2)<5(ri - r2)(5a/3 

- + eA(ri)) + /iscr ■ B(ri)]yG^/3(ri,T’2;Ti - T 2 ) - y] At,^(ri)F^/3(ri,r 2 ;n - T 2 ), 

7 7 

- T 2 ) = 

+eA(ri)) +Mscr-B(ri)]a^F^/ 3 (ri,r 2 ;Ti - T2) - E Aq/^y (t’i)G 7 / 3 (t’i,T 2 ;ti - r 2 ), 

7 7 

-9ri-F’„;3(ri,r2;ri - T 2 ) = 

- H-eA(ri)) +/is<7 • B{ri)]l^F^i3{ri,r2;Ti -T 2 ) -'^Al^{ri)G^i3{ri,r2;Ti -T 2 ) (B6) 

7 7 


where T,- denotes the imaginary time ordering operation, 
and 

c„(r) = 

c<^(r) = (B5) 

Here, /x is the chemical potential, and Af = J2k a 
is the particle number. By taking derivatives of the 
Gor’kov Green’s functions with respect to the imaginary 
time, the following left- and right-sided Gor’kov equa¬ 
tions are obtained: 


and 


dr2Gap{ri,r2; n - T 2 ) =(5^(ri - r2)5{Ti - T2)Sap 

+ '^Ga^{ri,r2-,Ti - T2)[C(-*V2 +eA(r2)) + HsCr • B(r2)]^/3 - E'^ 2 ;n - r2)A;^^(r2), 

7 7 

dr2Gap{ri,r2;Ti - T 2 ) =S^{ri - r2)S{Ti - T2)5ap 

- ^(5a7(r-i,r2;ri -T2)[C(V2 + eA(r2)) +/iscr • B(r2)]4 “ E (ri,r2;Ti - r2)A^/3{r2), 

7 7 

5x2-Fa/?( t*!, f' 2 ; n - T 2 ) = 

-'^Fa.y(ri,r2;Ti - T2)lC(iV2 + eA{r2)) + iiscr ■ B(r2)]4 - E '^ 2 ; n - T2)Ajp{r2), 

7 7 

dr2Fap{ri,r2;Ti - T 2 ) = 

E (ri,r2;ri - r2)[C(-*V2 + eA{r2)) + HsCr ■ B(r2)]^/3 - E G {ri,r2-,Ti-T2)Al^l^{r2), (B7) 

7 7 


where — e is the electronic charge, A is the vector poten- and 
tial associated with 

^{k) = {sk-l^)(To+Cgk-CT. (B9) 


Aap{r) = {-ia2)apA{r) 


(B8) Hereafter, we define the operation of V to an arbitrary 
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function h{r) from the right side as 
h(r)V = -Vh(r). 

The Wigner representation of the Green’s functions is 

r , 

/ ^3^/g-ifc.r / dre^^-^Xc,^{r + r'/2,r-r'/2;T) 

(Bll) 

for X = F. By Fourier-transforming Eqs. dEH) 


and (lB7l) and neglecting the higher-order terms with re¬ 
spect to l/k^^x and l/k^^z {^x and are the coherent 
(BIO) lengths in the x and 2 ; directions), the left- and right-sided 
Gor’kov equations in the Wigner representation become 


G-^G = GG-^ = i 

in the matrix form, where 

\ - f G{r-,k,ujn) F{r-,k,u)n) \ 
-G{r-,k,ajn)) 


(B12) 


(B13) 


G ^(r;fe,w„) = 

GunCro - [^(fc) 


\-v{k)-11/2 +HsCr ■ B{r)] ia 2 A{r) 

-ia 2 A*{r) -iw„(7o - + v'^{-k) ■ Tl*/2 + /igcr"’’ • B{r) 


(B14) 


and i is the 4x4 identity matrix. Here, 

n = -iV + 2eA, 
n* = iV + 2eA, 


X = 


X-f-f X-fi 

Wt Wt 


iX = G,G,F,F) 


and 


v{k) = VfcCW- 


(B15) 

(B16) 

(B17) 


The subtraction of the two equations in Eq. (jB12p leads 
to the following equation 


[G-\G]=0. 


(B18) 


To obtain Gor’kov equations in a more useful form, the 
transformation with the matrix 


0 u{-ky 


(B19) 


is considered, namely, the Green’s functions and its in¬ 
verse operator in the transformed new representation are 
defined as 


G' = tJ{k)GU\k), 


and 

G'-i = U{k)G-X\k), 
and further, we put 
/ Gr 


G' = 


G 21 


Gi2 

Fi 

7^12 \ 

G2 

F 21 

F 2 

— Fi2 

-Gi 

-Gi2 

-F2 

-G 21 

-G 2 / 


(B21) 


(B22) 


Suppose here that the length scales on any inhomogene¬ 
ity are sufficiently longer than k^^. Then, Xi and X 2 
are interpreted as the intraband Green’s functions of the 
bands 1 and 2 at r, while X 12 and X 21 are interpreted as 
the interband ones {X = G^G^F^F). Here, we neglect 
the off-diagonal elements by assuming that Tc^ the critical 
temperature at zero field, and the magnitude of the 
Zeeman energy, are much smaller than C, the strength of 
the ASOC (see the main text). Then, Eq. (jB18[) becomes 


[G^ Ga] — 0, 


where 


Ga = 


Ga Fa 
-Fa -Ga 


(B23) 


(B24) 


(B20) and 


Ga = 


- [Ca + Va * 11/2 -h (-1)"+Vs^fc ' B{r)] 
-<A*(r) 


WaA{r) 

-i^n - [ia - Va ' n*/2 - (-1)"+^^^^ * B{r) 


(B25) 
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Im 



FIG. 12: Two paths of the complex integration f. 


Here, 


Ca — ^ak ^5 


(B26) 




(B27) 


and 


(z zj =tew ^;;w) 

In this representation, we define the quasiclassical 
Green’s functions on each FS as follows: 


Appendix C: Conditions on Qa 

Because of the subtraction which leads to Eq. (jBlSp . 
some information is lost. Actually, we cannot determine 
the quasiclassical Green’s functions uniquely based only 
on Eq. (|B34p . This information, however, can be recov¬ 
ered with the following two conditions: 

5a = t sgnRe^a =-sgnw„, (Cl) 

where 1 is the 2x2 identity matrix. In this section, we 
derive these conditions based on some assumptions. 

According to the Eilenberger equation (jB34p . we have 
the relations 


'^Fa 


■ V(fi(a - 5a) = 0, (C2) 

-/Ta|=0, (C3) 


where Vpa = 

In the spatially homogeneous case, the solutions of the 
Gor’kov equations in Eqs. (jB6p or in Eqs. (IB7p lead to the 
quasiclassical Green’s functions 


fa 


Qa — Qa — 

WgA 

W'^n + teaA|2 


LUy] 


yjujl + Iw^aAP’ 

{WgA)* 


,fa = 




,A|2’ 


(C4) 

(C5) 


which obey 


5a(r;fc,w„) = j 

f —Gair-,k,UJn), 

TTl 

(B29) 

9a{r;k,u}n) = j 

1 —Ga{r-,k,UJn), 

TTl 

(B30) 

fa{r;k,uJn) = < 

fd^a ^ 

p — rFa{r-,k,UJn), 

1 TTl 

(B31) 

fa{r-,k,Un) = 1 

i^Fair-,k,LOn), 

1 TTl 

(B32) 


and their matrix form as 

S. = (4 4) ■ (B33) 

Here, the complex integration f represents the average 
of the two contour integrals along the paths 1 and 2 il¬ 
lustrated in Eig.[T2l The integration of Eq. (jB23p with 
respect to leads to the Eilenberger equation: 


5a = i- (C6) 

These are the physical solution of the Eilenberger equa¬ 
tion (IB341) . 

In the cases with spatial inhomogeneity, it is assumed 
that the system can be smoothly transformed to the ho¬ 
mogeneous state towards infinity so that the normaliza¬ 
tion condition (jC6p is still valid according to Eqs. (jC2p 
and (jC3p . Therefore, generally, 

9a=9a = -\/l + fala (C7) 

is also valid. Strictly speaking, the branch of the root in 
Eq. (jC7j) cannot be determined based only on Eq. (jC6j) . 
However, we note that the inequality 

I/a/a I < 1 (C8) 


G 


-1 I 


\k = kFa 


^Qa 


= 0 , 


(B34) 


where k-pa is the Eermi wave vector of the ath band, and 
we also use the fact that every Green’s function has a 
sharp peak with the width |A| around = 0- 


is satisfied in the homogeneous case, and that the aver¬ 
aged magnitude of the superconducting energy gap |A| 
becomes smaller due to a spatial inhomogeneity of A. 
Thus, the inequality (jC8j) should remain valid in vortex 
states so that the same branch as in the homogeneous 
case may be chosen. 
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Appendix D: Symmetry Relations of Quasiclassical 
Green’s Functions 


In this section, we derive some symmetry relations con¬ 
necting one quasiclassical Green’s function with another. 

In the Wigner representation, the Gor’kov Green’s 
functions have the symmetry relations 


-^a/3 ^n) 

Fap{r-,k,UJn) 
Fap{r-,k,UJn) 
Gocjiix 1 ^n) 

Gap{r-,k,Un) 
Gap{r-,k,Un) 


k^ ^n)i 
-Fpa{r-, -k, -UJn), 
Fpa{r;k,-ujn)*, 
-Gpa{r-, -k, -UJn), 
Gi3a{r-,k,-L0n)*, 
Gi3a{r-,k,-L0n)* 


(Dl) 



following from their definition. In other words, with the 
use of the transformation (IB2QII . we have 


Fa{r-,k,Un) 

Fa{r-,k,Un) 

Fa{r-,k,Un) 

Ga{r-,k,Un) 

Ga{r-,k,Un) 

Ga{r-,k,Un) 


-Fa{r-, -k, -Un), 
-Fa{r-, -k, -Un), 
Fair; k,-UJn)*, 
-Gair;-k,-Un), 
Gair;k,-uJn)*, 
Gair;k,-uJn)*. 


(D2) 


Integrating these equations with respect to and using 
the fact that ga = Qa lead to 

fair; k, UJn) = - fair;-k,-UJn), 
fair; k, UJn) = - fair;-k,-UJn), 
fai‘^',k,UJn) — fai‘^',k, UJn) , 

gair;k,uJn) = -Qair;-k,-Un), 
gair;k,uJn) = -gair;k,-ujn)* , (D3) 

from which useful relations 

fLUn) — 5 

ga{r;k,uJn) = ga{r;-k,ujny (D4) 


are obtained. 

We often use these relations in this paper when sum¬ 
mands in k or uJn summations include the quasiclassical 
Green’s functions. 


Appendix E: Approximation on Fermi Velocity 


FIG. 13: Relations among the three vectors in Eq. (lE4l) . FSo 
denotes the Fermi surface of the bare band, while FSa’s (a = 
1, 2) express the two Fermi surfaces split by the ASOC of 
Rashba type. In the figure, FSa is represented by that of 
band 2. 


and 




k± = k±/ |fe_L|. 

(E3) 

Here, we put 




fepo = kpo + Ska, 

(E4) 


where k^a and k^o are the Fermi wave vectors 

on the Fermi surfaces split by the ASOC and that on the 
bare band, respectively (see Fig. nil. Keeping 


kp \EfJ 

(E5) 

in mind, we get 


Vp ■ Ska = -i-l)‘^+\ 

(E6) 

from 'Hsingie: whcrc t;F = '^0 • Hereafter, the terms 

ofO((C/EF)2), 0 ((J/F;f)^) andO(JC/E|) are neglected. 
In this approximation. 

Ska = -(-l)“+^TOTfe_L, 
Kf 

(E7) 


Substituting this expression to Eq. (jE4p and using 
Eq.jEl lead to 


VFa = '^F- (E8) 


The velocity in the band a is given, following its defi¬ 
nition, by 


Va 




where 


(El) 


(E2) 


Appendix F: Derivation of Free Energy 

The free energy measure from that in the normal 
phase, i.e., 

F = -Tlntre”'^'^ +Tlntre“^’^l^=“ 


1^0 = VfeEfe, 


(FI) 
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(y5 = 1/T), is used in the text to determine the vortex 
lattice structure at each field and temperature. Since ob¬ 
taining a tractable expression of the free energy directly 
from the above expression is not easy, it is first rewrit¬ 
ten by following the procedure based on the variational 
principle [H used by Eilenberger. 

For this purpose, the gap equation and the expression 
of the electric current density are needed. According to 
the definition of the order parameter (lB3l) and the rela¬ 
tion (lD4l) . the gap equation is 

-A + 2nT ^(K/a)FS = 0, (F2) 

^ 0<UJri<^c,CL 

where Na is the normal DOS on each FS, 



is the average over each FS for an arbitrary function h(/c), 
and uJc is the frequency cutoff introduced to prevent the 
divergence of the summation. To treat uJc implicitly, we 
transform this equation to 


where 

js = + eA)cp{r))^^, (F7) 

a,l3 


3m = -mW X {ci{r)(Tapcp{r)) . (F8) 

a,(3 

The relative current density components measured from 
their normal counterparts are expressed in terms of the 
quasiclassical Green’s functions as follows: 

Ajs =js - js'lA =0 

k,UJn,Oi,l3 

= - ie2TTT Y {VaiOa + 1))fs > (^9) 

UJn>0,a 

^3m =jM - 

=V X AMp^ra, (Fio) 


77111 ^A + 27rr Y 

^ UJn>0,a 



iw*afa 


+ 



= 0 


(F4) 

by introducing the mean-field transition temperature Tc 
at zero field through the well-known relation 


- - 2nTN Y 

g ^ ^ 


1 T 

— = A^ln—, 


(F5) 


T 

Z\ATpara — ^ ^ ^ ^ ^ 

k,ujn ,a,g 

= -ilis2T:T Y -^o((-l)“’^^9fe(5a + l)>Fs- 

UJn>0,a 

(Fll) 


Here, 


Here, N is the average of the normal DOS on the two 
FSs. The current density is obtained from the relation 

3 = - (F6) 

_I 


AGpa{r;k,U}n) = Gpa{r-,k,LOn) - Gpa{r-,k,LOn)\^^Q 

(F12) 

and we have used the fact that 5'a|A=o = for cjn > 0. 
Then, we define the expression 


n[A, A, /, f]=- J d^r {Ajs ■ A + AMp^^a • B) 

Na 


/ 


+ / d^r 


N\Af\n- + 27TT Y 


UJn>0 


iA*W*Ja + iWaAfa + 


\WaA\ 


{da + 1) f 


2uJr, 


-i;f • Vlni^ 

2 fa 


FS 

(F13) 


as the functional from which the Eilenberger equations 
m and (HU), the gap equation (lF4l) . and the difference 
of the current Ajs + AJm follow after variations with 
respect to fa and /a. A*, and A, respectively. Next, 


by replacing fa and fa with the solutions /^[A, A] and 
/a [A, A] of the Eilenberger equation under given A and 
A, Eq. (|F13p is rewritten in the form 
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n = n [a.a./ia.ai./ia, A]] = / <i> 


iv|Api„i: + 2.r V + 

Tc 2 \ \ OJn 5a 1 / / FS 


r 


ilM 


obeying the conditions 

(5f2 _ sn 

5A~ 5A 


sn _ 6^ 

SA~ SA 


_ 6F 

/<,=/a[A.A] “ JA’ 
/„=/„[A,A] 


- ^ 
/a=/a[A,A] “ 5 A' 
/„=/„[ A, A] 


(F15) 


(F16) 


fi[A = 0,A = 0 ] = F[A = 0,A = 0 ] =0. (F17) 

Thus, 

F = n, 

which coincides with Eq. (j27j). 

Appendix G: Calculation of 


(F18) 


Here, the relation (j22]) is derived. 

For ujn > 0, 

[2 {w„+i(-l)“+Vs9fe ■-B} + inp • n] ^ 

= /“dpe-2K+F-i)“+V39fe-B}pg-i„p.np^ (Ql) 

The operators 

a =^(7“'/'nQ,, - n'/'nQ,,), (G2) 

=^(7“'/"nQ,, + n'/'nQ,,), (gs) 

which fulfill the relation [a, a’*'] = 1, are the annihi¬ 
lation and creation operators of the LLs ([T9]), where 
7 = e./e., Q = ^SNQo {SN = {N2 - Ni)/{Ni N 2 )), 
vh = ll\/2eH^ and 


IIq — II + Q. 


(G4) 


With the identity in the case where 

[A, [A,H]] = [H,[A,H]]=0, 

^-ivF-TlQp _ p‘^/2^-is*pa^ ^-ispa ^ 

Thus, with the definition of the Ath LL '0Ar(^) (see 
Eq. (fT^ in the main text), 

p), (G6) 


where [ • ]uc is the average over the unit cell and the 
relation 


[V’mV’JvIuC = ^M,N 


is used. Here, 


s = 




(G7) 


(G8) 


^. (.)"-(-2r-' 


1=0 


(M-/)!(A-/)!/! 


(G9) 


From Eqs. (jGip and (jG6p . 

[2 {un + i{-lT^^fisgk • • n]“^A 

— 2|u;ri,+h~l)“"^^Ms^fc--B}PgivF-Qpg—iVF-HQp^ 

= (GlO) 

with the definition of the matrix 

poo 

Jo 

X (Gil) 

which leads to Eq. ( 


Appendix H: Derivation of LDOS 

Here, we define the retarded Green’s function as usual: 

G^^(ri,r2;ti - ^ 2 ) 

= - t2) (||ca(ri,ti),cyr2,t2)|)^^ , (HI) 

where .b| = AB + BA is the anticommutator of 
arbitrary operators A and B, 

'l (t>0) 


m = 


(H2) 


0 {t < 0) 

is the step function, and 

X{t) = pm-P-M) (JJ3) 

is the Heisenberg representation of any operator X. Its 
Wigner representation is defined as 


yoR 

^a/3 


B + OO 


iEt 


{r-,k,E)= f ( dte 

J J —00 

xG%{r + r’/2,r-r’/2;t). (H4) 
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As is well known, this can be obtained from Gaf 3 {r; k^ujn) Finally, after performing the unitary transformation 
by the analytic continuation E ^ irj is an in- Eq. (jB2Qp in the above expression, Eq. (|32)) is obtained, 

finitesimal positive number). 

The LDOS N{r]E) can be defined by using this func¬ 
tion. 


N{r-E) = --^Y.^mGl^{r-,k,E). (H5) 

TT V ^^ 

k,a 
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